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We evolve a linearized (Teukolsky) solution of the Einstein equations with a non-linear Einstein 
solver. Using this testbed, we are able to show that such gravitational waves, defined by the Weyl 
scalars in the Newman-Penrose formalism, propagate faithfully across mesh refinement boundaries, 
and use, for the first time to our knowledge, a novel algorithm due to Misner to compute spherical 
harmonic components of our waveforms. We show that the algorithm performs extremely well, even 
when the extraction sphere intersects refinement boundaries. 


I. INTRODUCTION 

Computing candidate waveforms for gravitational 
wave observatories remains one of the premiere goals 
of numerical relativity. In order to achieve this goal, 
one needs to sensibly define radiation, develop numer- 
ical codes and algorithms to solve the non-linear Ein- 
stein equations, and handle a variety of technical prob- 
lems such as modeling singular black hole spacetimes and 
allocating available computer resources in problems with 
multiple dynamical length scales. 

In this brief report, we define a simplified problem in 
which to study the propagation of gravitational waves 
and to study a candidate algorithm due to Misner [1] 
for compute the spherical harmonic components of radi- 
ation fields. Consistent with our recent efforts to intro- 
duce mesh refinement techniques into numerical relativ- 
ity, to help deal with the problem of efficiently allocating 
computer resources in problems with multiple dynami- 
cal length scales, we use this test problem to study how 
faithfully gravitational waves propagate across mesh re- 
finement boundaries, and to validate the Misner algo- 
rithm in the case where the spheres on which the spher- 
ical harmonics are computed cross through refinement 
boundaries. 

We are motivated to consider the case where extrac- 
tion spheres cross refinement boundaries in part because 
in our particular fixed mesh refinement scheme (provided 
by the PARAMESH package [2]) we find that there are fre- 
quently few or no radii for extraction spheres inside sin- 
gle mesh refinement regions when we place refinement 
boundaries at locations that most efficiently use avail- 
able computer memory. More importantly, however, we 
look forward to using adaptive mesh refinement, in which 
we anticipate having refinement boundaries dynamically 
pass through extraction regions as high resolution grids 
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track wavefronts. In this case, it will be essential that 
our spherical harmonic decomposition algorithm work in 
the presence of refinement boundaries. 

Finally, we note that setting down in detail a simple 
but useful test case can play an important role in making 
comparisons between different codes and different algo- 
rithms, much in the spirit of Ref. [3]. This is especially 
useful for the spherical harmonic problem since, at this 
point, few if any other numerical relativity groups have 
begun implementing this technology. We hope that this 
work can serve as a benchmark for future efforts. 


II. METHODOLOGY 

Computing gravitational waves from numerical simu- 
lations of the Einstein equations requires combining a 
variety of nearly mathematical formalisms and numer- 
ical techniques into a single code. In Section II A, we 
briefly discuss our particular Einstein solver. We are de- 
liberately brief on this point because, for our present pur- 
poses, we find it useful to abstract away from the details 
of how the simulation is done, and consider the problem 
of analyzing data once we have it. The methods that we 
present, implement, and test here are not specific to any 
particular formalism or to our code. These more general 
issues are made explicit in the rest of the section where we 
define gravitational radiation via the Weyl scalars in the 
Newman-Penrose formalism in Section II B, and discuss 
the need and a method for computing spherical harmonic 
components of radiation data in Section II C. 


A. The Hahndol Code 

Our code is a fully three-dimensional, non-linear Ein- 
stein solver. We use a fairly standard formulation of the 
3+1 Einstein evolution equations known as BSSN [4, 5]; 
our particular implementation was described in detail in 
Ref. [6] . Because the formulation of the equations is not 
of primary interest here, and because the basics of the 
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BSSN system are widely know, we do not focus on these 
equations. 

The code uses second-order accurate finite differenc- 
ing to approximate spatial derivatives and the iterative 
Crank-Nicholson method [7] to integrate the evolution 
equations forward in time. In both our previous and cur- 
rent work, we have verified second order convergence in 
our results. 

Our code parallel code uses the PARAMESH libraries [2] 
to handle domain decomposition and inter-processor syn- 
chronization. In addition, these libraries enable us to use 
non-uniform grids to focus computational resources in 
specific areas of the computational domain. Although 
the libraries support, and we look forward to using, the 
ability to adaptively modify the grid structure during 
the course of a simulation (adaptive mesh refinement), 
we currently fix our grid structures in advance using a 
priori estimates of where to focus resolution (fixed mesh 
refinement). 

Previous studies, showed that the Hahndol code is able 
to propagate strong gravitational waves (defined in terms 
of metric components) across mesh refinement bound- 
aries [ 8 ], and that the same code can handle strong, 
dynamical potentials moving between refinement regions 
[ 6 ]. In this paper we again focus on wave propagation, 
this time using a more formal definition of graviational 
radiation, the Weyl scalars in the Newman- Penrose for- 
malism, and on analyzing such data in a meaningful way. 


B. Weyl Scalars 

We use the Newman-Penrose formalism to compute 
gravitational radiation quantities. In this formalism, one 
chooses a tetrad of four null basis vectors. These vectors 
are conventionally labled l a , n°, ra a , and m a . 1 Of these, 
l a and n a are real and point along “outgoing” and “in- 
going” directions, respectively. The vectors m a and m a 
span “angular” directions and are complex conjugates of 
each other. These vectors are always chosen to satisfy 
the orthogonality conditions 

l a m a = n a m a — 0. (1) 

In addition, we impose the normalization conditions 
/ a n a = —1 and m a fh a = + 1 , which are imposed by most, 
but not all, authors (cf. Ref. [9]). 

Given a tetrad, tensor quantities can be recast as sets 
of coordinate scalars by projecting tensor components 
onto the basis. Applying this procedure to the Weyl ten- 
sor C pqrs yields five complex scalars. This is particularly 


1 Note that, although our code uses a 3+1 decomposition of the 

Einstein equations, indices on Newman-Penrose quantities are 
always understood to run over four dimensional spacetime in- 
dices. 


useful for gravitational wave studies since one of the five, 

^4 = —C PQrs n p fh q n r rh s , (2) 

represents outgoing graviational waves. Indeed, the pri- 
mary goals of this paper are to construct ^4 in some sam- 
ple, numerically evolved spactimes, and to demonstrate 
that such computations are accurate characterizations of 
the wave content of the spacetime. 

In vacuum, which is the only case that we consider 
here, the Weyl tensor is numerically equal to the Rie- 
mann tensor. 

As a technical point, it is important to note that, al- 
though our code uses a 3+1 decomposition of the Einstein 
equations, Eq. 2 is written in terms of four dimensional, 
spacetime quantities. Since, in code, we have only 3+1 
quantities on a single time slice from which to compute 
the four dimensional Riemann tensor, we follow Ref. [10] 
and write 

Rijki = Rijki + 2Ki[ k K(]j (3a) 

{4) Ro jkt = -2(K m] +T p j[k K l]p ) (3b) 

(4) i?0j0i = Rji ~ K jp K p t + KKji (3c) 

which express four dimension quantites on the left hand 
sides in terms of 3+1 quantities on the right . 2 Here 
Kij is the extrinsic curvature tensor for a spatial slice 
as imbedded in the full spacetime manifold, and K is its 
trace. The symbols r? fc , Rijki t and Rji stand for, re- 
spectively, the connection coefficients, Riemann tensor, 
and Ricci tensor associated with the three dimensional 
spatial metric on the slice. 


C. Spherical harmonic decomposition 

In our simulations, we would like to be able to ex- 
tract the spherical harmonic components of gravitational 
waves. We find this very valuable when analyzing the 
data because 

1 . Numerical noise tends to have higher angular fre- 
quency than genuine wave signals, and is therefore 
filtered by the decomposition process. 

2 . A priori knowledge about symmetries in the data or 
dominant modes associate with physical processes 
allow important checks on the plausibility of nu- 
merical solutions (especially when exact solutions 
are not available). 


2 There are two errors in Ref. [10] associated with what we call 
Eq. 3. First there is a factor of two difference between Eq. 3b, 
which is correct, and the corresponding equation in Ref. [10]. 
Second, Eq. 3 in this paper properly reflects the fact that the 
left hand sides refer to the four dimensional Riemann tensor. 
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3. Analysis of spherical harmonic components will be 
important in quantifying the amount of energy and 
momentum carried by gravitational radiation. 

In practice, however, we face a technical problem in 
computing spherical harmonic components of our data. 
While our data is stored on a cubic grid, the spherical 
harmonic components 

r) = j s F; m (0, r, 9, (p)dQ (4) 

of a general (complex valued) function $ in a spherical 
harmonic basis with spin- weight s are defined by an in- 
tegral over a sphere. Computing the integral in Eq. 4 
requires some type of interpolation from grid points to 
points on the sphere of radius r. 

One could solve this problem by the straightforward 
method of (1) interpolating the grid function $ to points 
on the sphere, and then (2) using some numerical ap- 
proximation to the integral in Eq. 4. This process of 
iterpolation followed by integration would then need to 
be performed at each time, at each radius, and for every 
function for which the spherical harmonic components 
are desired. 

We adopt a different algorithm due to Misner. Follow- 
ing Ref. [1], we smear the surface integral in Eq. 4 into 
a volume integral over a shell of half-thickness A, and 
create an orthonormal basis for functions on this shell by 
combining the spherical harmonics in the angular direc- 
tions with the first N Legendre polynomials in the radial 
direction. This approach, which we describe in more de- 
tail in Appendix A, has the advantage that, given the 
parameters A and AT, one needs only compute a rela- 
tively small number of weights for computing the volume 
integral, and that these weight depend only on the grid 
structure and the radius of extraction. This means that 
the weights can be computed once at the beginning of 
a simulation and stored. Further interpolations are not 
need, and the weights are valid for all functions. 

The question of how to choose the input parameters 
A and N was addressed in Refs. [11, 12]. Under the 
assumption (motivated by the analysis in the references) 
that A is chosen proportional to the grid spacing, the 
parameter N controls the order of convergence of the 
decomposition algorithm, while the parameter A controls 
the size of the error at that order. This analysis leads to 
the 

Rule of Thumb: Choose N just large 

enough to ensure that the error term pro- 
portional to A is not of leading order in grid 
spacing. Choose A just large enough to safely 
resolve Pn on the shell. 

For second order accurate codes, like ours, with grid 
spacing h , Refs. [11, 12] suggest choosing N = 2 and 
A = 3/i/4 based on this rule and emprical experiments. 
In cases where the shell passes through multiple re- 
finement regions, the grid spacing of the coarsest grid 


through which the shell passes should be chosen to en- 
sure that Pn is resolved on both sides of the interface. 
We follow these suggestions for all work presented here. 


III. TEUKOLSKY WAVES 

We apply our techniques to two sample problems. The 
first is a linearized solution to the Einstein equations due 
to Teukolsky [13]. We choose the Teukolsky wave as our 
first test problem because (1) it has an exact solution, 
so we can compute our numerical errors exactly, (2) it is 
linear, so we can extract our waves at small radii. The 
second point saves time and resources in the early stages 
of testing. Waiting for waves generated by non-linear 
sources to propage to the wave zone is costly, and cou- 
ples problems with evolving the sources to problems with 
extracting the waves. 

We tackle a non-linear problem in our second test case, 
a head-on collision of equal mass black holes, in Sec- 
tion IV. 


A. Analytic Preliminaries 

Our first problem studies a Teukolsky wave spacetime 
[13]. This is a solution to the linearized Einstein equa- 
tions, which represents a weak gravitational wave prop- 
agating through space. The linear nature of the initial 
data, makes this an excellent first test problem for two 
reasons. First, the there is an analytic solution for all 
times, allowing an exact calculation of numerical errors 
for convergence and accuracy studies. Second, because 
the wave is linear in the initial data, we are able to ex- 
tract our waveforms at small radii and short evolution 
times. This second fact is extremely useful for debugging 
algorithms, and also allows higher resolutions for a fixed 
problem size since the boundaries of the computational 
domain need not be as large as it would in problems (like 
the head-on collision described in Section IV) in which 
the waves are generated by non-linear sources and must 
propagate to the wave zone before being extracted. 

Although the Teukolsky solution is well-known, we 
summarize it here for completeness and to establish no- 
tation. We use the pure l = 2, m = 0 (in a spin-weight 
—2 basis) solution. 

The spacetime metric, given in full in Appendix B, is 
written in terms of functions 
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and in terms of a free generating function F = F(t - r). 
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The notation 


p( n ) _ 


’(j n F(x)' 

dx n 


x=t—r 


( 6 ) 


denotes various derivatives. Taking F as a function of 
t — r corresponds to outgoing waves. To generate an 
ingoing solution, change the argument of F to t + r, and 
change, in Eq. 5, the sign in front of all of the terms with 
odd numbers of derivatives. 3 

Furthermore, we follow Choi et al. [8] in choosing 


For this relatively simple spacetime, one can compute 
the value of the Newman- Penrose quantity 


*4 = 


sin 2 9 
16 


d 2 C d 2 A (d 3 B d 3 A\ 

~ 12 ~W +6 W + r I 3 — - + ^> 


dt 3 dt 3 J J 


( 9 ) 

which represents outgoing gravitational radiation. Not- 
ing that 


-2>20(M) = \/ 3^ Sill26f 
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(10a) 


F(x) = g e-W 


7 ) 



Y 00 (6,4>) ~ ^Y 2O (0,4>) 


10b) 


as the exact form of the generating function, where the 
free parameters A and A represent the amplitude and 
the width of the wave respectively. The natural length 
unit in the problem is A, and we consistently choose A = 
2 x 10“ 6 A. Moreover we take an equal superposition of 
an ingoing wave and an outgoing wave, both centered at 
the origin, for the initial data. This particular choice has 
a moment of time symmetry that allows us to set the 
extrinsic curvature tensor to zero. 

Choosing F of the form in Eq. 7 gives a waveform with 
oscillations but of essentially compact support. This is 
ideal for testing codes with boundaries because it makes 
clear when the wave passes through those boundaries, 
and allows one to easily detect any reflections that occur 
due to poor interface conditions. 

For our 3+1 code, we use this analytic solution to gen- 
erate the intial data, and evolve forward in time using 
the gauge conditions lapse a = 1 and shift ft 1 = 0. 

In this work, we use the Kinnersley tetrad [14] 

l — ~(r 2 -Fa 2 , A,0,a) (8a) 

n = +(r 2 + a 2 , —A, 0, a) (8b) 

m = — -= - (iasin0, 0, l,i csc0) (8c) 

V2(r 4- ia cos 6) 

for the linear wave problem. 4 Eqs. 8 give the tetrad 
in terms of the Boyer-Lindquist coordinates for a Kerr 
black hole. In these equations, A = r 2 — 2 Mr + a 2 , 
E = r 2 + a 2 cos 2 0, M is the mass of the Kerr black hole, 
and a is its spin. Because Boyer-Lindquist coordinates 
are compatible with the coordinate system used for the 
Teukolsky metric (Eq. Bl) when M = a = 0, we find this 
coordinate expression for the tetrad completely accept- 
able for this problem. 


3 Note that Ref. [8] improperly describes the construction of in- 
going waves. The description here, which matches the original 
reference, Ref. [13], is correct. 

4 The first component of m a should have the sin# term in the 

numerator as in Eq. 8c; the corresponding equation in Ref. [10] 

is incorrect. 


it is clear that, as claimed, Eq. 9 is a pure l = 2, m = 0 
mode. Note also that Eq. 10b provides a way to com- 
pute a spin-weight —2 spherical harmonic component 
from the usual (spin-weight 0) spherical harmonic com- 
ponents. This is essential for us because, although we 
would like spin- weight —2 results, our current implemen- 
tation of the Misner algorithm computes only spin- weight 
0 components. 

B. Numerical Results 

We evolved initial data specified by the Teukolsky so- 
lution described in Section III A on a domain with outer 
boundaries at 32A and with a “box-in-box” refinement 
scheme. The boundaries were at 2A, 4A, 8A, and 16A; 
neighboring regions differed in resolution by a factor of 
two, with the finest regions surrounded by coarser re- 
gions. In order to compute convergence factors, we ran 
the code at innermost resolutions A/12, A/24, and A/48. 

We then attempted to compute the spherical harmonic 
components of that wave at five distinct “detectors” lo- 
cated at r = 3A, 4A, 5A, 6A, and 7A. To help visualize 
how the refinement boundaries and the extraction radii 
are related geometrically, we find it useful to draw ex- 
traction maps , as shown in Fig. 1. The extraction maps 
make it clear that the shells used to compute the spher- 
ical harmonic components generically pass through mul- 
tiple refinement regions, especially since in three dimen- 
sions the corners of the cubic refinement regions tend to 
pass through the spherical extraction shells. The extrac- 
tion maps shown in Fig. 1 correspond to the innermost 
(r = 3A) extraction radius in the coarsest run, and the 
r = 6A extraction radius in the medium resolution run. 
(Extraction maps for nearly all extraction radii and two 
resolutions appear in Ref. [11].) Note that only the in- 
nermost refinement regions are shown in the maps. In 
each case, the entire map was surrounded by additional 
(coarser) refinement regions. 

Since we know the analytic solution for this test prob- 
lem, we are able to make a detailed convergence study 
with three resolution levels. We were able to verify 
that not only the waveforms themselves, the extracted 
spherical harmonic components converge at second order 
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FIG. 1: Selected extraction maps. Each map is labeled with 
a triple of numbers (TV, R, A), which indicate the number of 
(cell-centered) grid points across one coordinate direction in 
the coarsest region, the extraction radius (in units of A), and 
the half-thickness of the shell. Note that the shells generically 
pass through multiple refinement regions, especially since, in 
three dimensions, the corners of the cubic refinement regions 
tend to poke through the spherical extraction shells. The 
three circles drawn on each graph show the extraction sphere 
and the edges of the finite-thickness shell around the sphere 
used in the Misner algorithm. 


(as predicted), even with the extraction spheres passing 
through the refinement boundaries, and that the solution 
is highly accurate. We also find that the results are ex- 
tremely consistent between the various extraction radii. 
Note in particular Fig. 2, which shows, in Panel A, a sin- 
gle wave as computed at distinct radii. The amplitude, 
to leading order, scales like 1/r and maintains its shape 
as it moves outward, consistent with the analytic solu-. 
tion. This is seen more clearly in Panel B in which the 
waveforms have been scaled-up by r and shifted to a sin- 
gle extraction radius. Up to higher order terms in 1/r, 
these waveforms should match. As seen in the figure, the 
agreement is excellent. 



Scaled and ahifted wave axiraction for TaukoUky initial data (a * 2e-6 A)on FMR jrkJ 



FIG. 2: The Z = 2, m = 0 component of the Teukolsky wave, 
as computed at five different radii (highest resolution). The 
wave form is preserved up to the leading order 1/r scaling. 
Panel A shows the raw data. Panel B shows the same data, 
scaled-up by r and shifted in time to t — 3A. 


In evaluating the effects of the refinement boundaries 
on the spherical harmonic decomposition algorithm, one 
should bear in mind that no two of the extraction radii 
have the same geometric relationship with the refinement 
boundaries (cf. Fig. 1), and yet they nonetheless generate 
perfectly consistent results. 


IV. BLACK HOLE HEAD-ON COLLISION 

We also tested our techniques on a non-linear prob- 
lem, which is more closely related to the astrophysical 
problems of interest to us. Our specific choice of the 
equal-mass head-on collision of two black holes tests our 
techniques on a fully non-linear problem, but also, con- 
sistent with our goal of testing our techniques, allows a 
variety of symmetries that we can use as checks on our 
numerical solutions. 


A. Analytic Preliminaries 

In our head-on collisions, we place two black holes, 
each of mass M, on the 2 axis at a coordinate distance of 
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1.1515M from the origin. The initial data is generated 
by the puncture prescription of Ref. [15], which is a gen- 
eralization of Brill-Lindquist initial data [16]. We then 
use the excision algorithm of Ref. [17] to avoid evolving 
the portions of our grid interior to the (seperate) event 
horizons of our initial data. We neither allow our exci- 
sion regions to move nor switch to a single, larger excision 
region after a common horizon forms. Although one or 
both of these techniques will likely be desirable in fu- 
ture simulations, we find the simple approach sufficient 
to extract convergent waveforms. 

Once we generate our initial data, we evolve it in our 
3+1 code using the 1+log evolution equation 


da 

~di 


-2 aK 


(ii) 


for the lapse a, and the hyberbolic Gamma-driver evolu- 
tion equations 


dp 

dt 

dB i 

dt 



(12a) 

dr' , 

IT-** 

( 12 b) 


for the shift /3\ In these equations, K is the trace of the 
extrinsic curvature tensor for our slice, T l is the confor- 
mal connection function of the BSSN evolution system 
(cf. Ref. [5]), and B l is defined by Eq. 12a to keep the 
gauge evolution equations first order in time. The quan- 
tity 


N 

i>BL(r) = ^2 

i= 1 


1 + 


Mi 


2 | r 


(13) 


is the Brill-Lindquist factor used to generate puncture 
initial data for N black holes with masses Mi and posi- 
tions Vi. (In our case, N = 2 and M\ = M2 = M.) In our 
work, we choose the dissipation parameter 77 = 2.8/M, 
and set, at the initial time, a = 1 and ft 1 — 0. These 
gauge conditions were first studied in Ref. [18]. 

For the head-on collision, we do not use the Kinnersley 
tetrad. We have no way of knowing for the head-on col- 
lision that the coordinate expressions in Eq. 8 are appro- 
priate to our numerically evolved spacetime. (Indeed, for 
a general problem, one would assume that they are not 
appropriate.) For the head-on collision, we instead con- 
struct an orthonormalized tetrad from the numerically 
evolved spacetime using a Gramm-Schmidt procedure as 
described in Ref. [10]. 5 


In the head-on case, we do not have an analytic solu- 
tion for the Weyl scalars, so we can not directly compute 
the error in our solution. The symmetries of the problem, 
however, provide us with several analytic checks on our 
numerical results. Before discussion these symmetries 
specifically, however, it is worth making explicit a re- 
lated point: there are three independent symmetry axes 
in this problem. The first is the symmetry of the physical 
problem, i.e. the axis along which the two black holes col- 
lide. The second is axis with respect to which the tetrad 
is computed. (Note, for example, the Kinnersley tetrad 
defined by Eq. 8, in which the existence of a prefered axis 
is manifest in the coordinate expressions.) The third is 
the axis with respect to which we compute the spherical 
harmonics. Only if we align all of these axes of symmetry 
(conventionally along the z axis) will all of our symmetry 
checks be true. 6 

Assuming that, as in the problem described here, all 
three symmetry axes are aligned, the numerical solution 
should have the following properties: 

• Re{\k 4 } is axially symmetric and is symmetric un- 
der the transformation z — > — z up to the round-off 
level. 

• Im{^4} converges to zero. 

• Viewed in a spin-weight —2 basis, the dominent 
contribution to Re{^ 4 } comes from the / = 2, m = 
0 mode. 

• Viewed in a spin-weight 0 basis, trunctation level 
errors in Im{^ 4 } appear in only in modes with odd 
values of / > 5 and even values of m. 


B. Numerical Results 

For the head-on collision simulations, we again use a 
box-in-box mesh refinement scheme. We place our outer 
boundary at coordinate distance 128M, and place refine- 
ment boundaries at 4 M, 8M, 16M, 32 M, and 64M (six 
levels total). We evole only one octant of our space- 
time, using appropriate symmetry boundary conditions 
to mimic a full grid. Because of the octant symmetry 
conditions used, only one of the black holes appears in 
our evolved numerical grid. (The other is accounted for 
by the symmetry boundary conditions.) We ran sim- 
ulations at three resolutions with interior resolutions of 
M/16, M/24, and M/32 in order to perform convergence 
tests. 


5 We use the Teukolsky wave as a check on our implementation 

of the Gramm-Schmidt construction. In the limit of perturbed 
flatspace, the Gramm-Schmidt procedure recovers the Kinnersley 
tetrad up to a trivial rescaling of l a and n a . Comparisions of 
waveforms extracted from Teukolsky wave spacetimes with the 
two choices of tetrad match very closely when this rescaling is 
taken into account. 


In addition, because, for example, the inner product 
§ - 2 ^ 22 ( 0 , 0)oVh(0, <p)dCl ^ 0 for all l > 2, it is impossible to 
convert between the spin- weight —2 and spin- weight 0 bases in 
problems lacking axial symmetry, unless one computes spherical 
harmonic components at all values of l in the spin-weight 0 basis. 
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Extraction Map: (64, 20, 0.375) 



FIG. 3: Selected extraction maps for the head-on collsion. 
The maps accurately show the relative sizes and positions 
of refinement regions and the extraction shells used by the 
Misner algorithm for extractions at r = 20 M and r = 40 M 
in our lowest resolution simulation. 


We excise a cubic region, centered on the puncture. 
The cube has sides of length 0.23125M The size shape, 
and location of the excision region remain fixed during 
the course of the run. 

In these simulations, we extract our waveforms at r = 
20 M, 30 M, 40M, and 50M. Extraction maps for the 
20M and 40M extractions in Fig. 3. Since this is a non- 
linear problem, we are forced to extract a larger radii 
than in the Teukolsky wave simulations in order to ensure 
that we are extracting in the wave zone. 

With our data from the three resolutions, we were able 
to evaluate the convergence behavior of our waveform. 
Fig. 4 shows a three-point convergence plot for the l ~ 2, 
m — 0 (spin-weight —2) component of 4> 4 extracted from 
our simulations. In the three point convergence graphs, 

we show (^4 HlS ^ “ *i Med) ) and 4(4>4 We< ^ — $4*°^) /9, which 
should coincide, up to the effect of higher order error 
terms, for our second order accurate code, and our choice 
of relative grid spacings. The agreement is impressive. 


Convergence plot of 4/4,20 for head-on collision goes 
here. 

FIG. 4: A convegence plot of the l = 2, m = 0 (spin-weight 
—2) component of 4> 4 for the equal mass black hole head-on 
collision problem. 

Scaled and shifted waveforms for the head-on 
collision. 

FIG. 5: A comparision between the l — 2, m = 0 components 
of 4/ 4 as computed at distinct radii. The waveforms have been 
scaled-up by r and shifted to the intermost extraction radius, 
r — 20 M. Up to sub-leading terms in 1/r, these scaled and 
shifted waveforms should coincide. 

As with the Teukolsky wave, we should also be able to 
demonstrate agreement betwee the waveforms extracted 
at distinct radii, and indeed our numerical results are 
consistent with this expectation. Fig. 5 shows the l = 2, 
m = 0 component of ^4 as extracted at four radii. As in 
Panel B of Fig. 2, the components have by scaled-up by r 
and shifted to the innermost extraction location for direct 
comparision. Up to sub-leading terms in 1 /r, waveforms 
scaled and shifted in this way should coincide, as they do 
in the figure. 


V. DISCUSSION 

We have presented a useful test problem for study- 
ing the propagation of gravitational wave through mesh 
refinement boundaries. We find, consistent with earlier 
results, that our interface conditions faithfully propagate 
gravitational waves. More importantly, we show for the 
first time that the Misner algorithm for computing spher- 
ical harmonic components on a cubic grid is viable in live 
simulations and validate in an application the numerical 
analysis performed on this algorithm in earlier work. Our 
results show that the algorithm works in the presence of 
mesh refinement boundaries even when the extraction 
spheres pass through those boundaries. 

Because mesh refinement regions tend to be cubes, we 
find in our current simulations that this last point in 
crucial. There is frequently little or no room to locate 
finite-thickness shells such that they remain entirely in 
a single resolution region when we use box-in-box fixed 
mesh refinement. Moreover, we expect that, in the fu- 
ture, when we run simulations with adaptive mesh refine- 
ment, that high resolution regions will track wave fronts. 
In that case it would be almost inevitable that refinement 
boundaries would dynamically pass through the extrac- 
tion region. 

These results pave the way for further studies of the 
propagation of waves through mesh refinement bound- 
aries and for applying the Misner algorithm in more as- 
trophysically interesting situation in which the waves are 
generated by simulated, dynamical sources. We hope also 
that the simple test problem that we have presented here 
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can serve as a standard testbed for comparing such tech- and that 
niques in the future. 

p(r; R, A) = r~ 2 5(R - r) (A6) 
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APPENDIX A: SUMMARY OF THE MISNER 
METHOD 


This appendix provides a more detailed look at the 
Misner algorithm for computing spherical harmonic com- 
ponents of a function represented on a cubic grid. Ad- 
ditional details can be found in Misner’s original paper, 
Ref. [1]. A detailed discussion of the tr unction error as 
a function of the algorithm’s parameters can be found in 
Refs. [11, 12]. 

Recall that the problem is to compute the spherical 
harmonic components defined by Eq. 4, of a func- 
tion <£ that is known only on vertices of a cubic lattice. 
(In this appendix, we suppress the spin-weight index s.) 
In order to begin, two definitions are need. First define 
a family of radial functions 


Rn(r-,R, A) = r 


- 


2n + 1 
2A 


Pn 


r-R 


(Al) 


in terms of the usual Legendre polynomials P n . Here R 
and A are parameters that will be associated with the 
radius at which the spherical harmonic decomposition is 
desired and half of the thickness of a shell centered on 
that radius. From this, define 


Ynim (r, 0, 0) = Rn(r)Yim(0 , </>) (A2) 


is a delta function. (Compare Eq. A4 to Eq. 4.) 

On a finite grid T, the inner product Eq. A3 will have 
the form 


(f\a) = E f( x )g( x )w x (A7) 

where each point has some weight w x . This weight was 
given the form 


w x - 


0 |r - R\ > A + h/2 

h 3 jr - R\ < A - h/2 

(A + h/2 - \R — r\)h 2 otherwise 


(A8) 

by Misner, where h is the grid spacing. Only cases with 
A > h/2 are considered. This means, roughly, that 
points entirely within the shell S are weighted by their 
finite volume on the numerical grid, points entirely out- 
side of the shell S have zero weight, and points near the 
boundary are weighted according to the fraction of their 
volume inside S, 


With the numerical inner product Eq. A7, and letting 
capital Roman letters A = (nlm) represent index groups, 
the Ya are no longer orthonormal. Their inner product 


(Ya\Yb)=Gab = Qba (A9) 


forms a metric for functions on the shell. (Although a 
priori this matrix appears to be complex valued, it is ac- 
tually real-symmetric and sparse, cf. Refs. [11, 12]. For 
now it suffices to follow Misner in denoting it as generi- 
cally Hermitian.) The inverse to this metric G AB can be 
used to raise indices on functions defined on the sphere. 

Making use of this new metric, and with some further 
analysis, the approximation for the spherical harmonic 
coefficients 


which form a complete, orthonormal set with respect to 
the inner product 

(f\9) = / f{x)g(x)d 3 x (A3) 

Js 


R) = ^2 Rim(x; R)w x $(t, x) 

x€V 


follows with 


(A10) 


on the shell S — {(r, 0, <t>) | r £ [R - A, R -f A]}. Note 
also that, because the functions R n form a complete set 
in the radial direction, 

R)= J p(r ; R, A )Y lm (6, 4>)$(r, 9, cfr)d 3 x (A4) 

with 

oo 

p(r; R, A) = J2 Pn(R; R, A)Rn(r; R, A), (A5) 

n=0 


N 

Rim(r-,R) = J^Rn(R)Y nlm (r,9,<P) (All) 

71=0 

in terms of Y A = G ba Yb-> not Y^. 

Note that one need only store the combination Wx&im 
at points where w x ^ 0 in order to compute the spheri- 
cal harmonic components. This buries all the details of 
the computation in a relatively small table of numbers, 
and, since these numbers are not time dependent, this 
calculation need be done only once per simulation. 
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APPENDIX B: TEUKOLSKY WAVE SOLUTION for the l = 2, m = 0 case. 


The general form of the spacetime metric for the 
Teukolsky wave solution [13] is 


ds 2 = — dt 2 + (1 + Af rr )dr 2 

+ 21 Bf r $rdrd8 + 2 Bf r< pr sin ddrdcj) 

+ (i + Cf$$ + r2 ^0 2 0^1) 

+ 2(A — 2 C)fe<f>r 2 sin Oddd(p 
+ (l + Cf$ + Af ffj r 2 sin 2 6W</> 2 


given in terms of the functions Eq. 5. The angular func- 
tion are 


f rr = 2 — 3 sin 2 0 (B2a) 

/ r $ = — 3 sin# cos# (B2b) 

fr<j> = o (B2c) 

/$ = 3 sin 2 0 (B2d) 

= -1 (B2e) 

fo* = 0 (B2f) 

4? = (B2g) 

4 ? = 3 sin 2 9 — 1 (B2h) 


The Weyl scalar ^4 for this spacetime is computed 
from the definition Eq. 2 using the Kinnersley tetrad. 
Eq. 8, and noting that, of the twelve non-zero compo- 
nents of the Riemann tensor associated with the metric 
Eq. Bl, only 


Rtdte 

Rt<j>t<j> 

RtdrO 

Rr(pr<f> 

RrOrO 


3 2 . 2 n&C 1 2 d 2 A 

= 

3 2 ■ 4 Jd 2 c 


(B3a) 


1 2 • 2 n &A 

+ 2 r sm d lF 


(B3b) 

1 0 , „ / d 3 4 5 6 7 8 9 B 

d 3 A\ 

(B3c) 

-8 r Sl ” * l 3 

+ dt 3 ) 

- sin 2 ORtete 


(B3d) 

sin 2 e Rt<t,t4> 


(B3e) 

- sin 2 ORtere 


(B3f) 


contribute to the sum. 
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